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We theoretically investigate the behavior of dark solitons in a one-dimensional spin-orbit coupled 
atomic Fermi gas in harmonic traps, by solving self-consistently the Bogoliubov-de Gennes equations. 

The dark soliton - to be created by phase-imprinting in future experiments - is characterized by a 
real order parameter, which changes sign at a point node and hosts localized Andreev bound states 
near the node. By considering both cases of a single soliton and of multiple solitons, we find 
that the energy of these bound states decreases to zero, when the system is tuned to enter the 
topological superfluid phase by increasing an external Zeeman field. As a result, two Majorana 
fermions emerge in the vicinity of each soliton, in addition to the well-known Majorana fermions at 
the trap edges associated with the nontrivial topology of the superfluid. We propose that the soliton- 
induced Majorana fermions can be directly observed by using spatially-resolved radio-frequency 
spectroscopy or indirectly probed by measuring the density profile at the point node. For the latter, 
the deep minimum in the density profile will disappear due to the occupation of the soliton-induced 
zero-energy Majorana fermion modes. Our prediction could be tested in a resonantly-interacting 
spin-orbit coupled 40 K Fermi gas confined in a two-dimensional optical lattice. 

PACS numbers: 03.75.Ss, 71.10.Pm, 03.65.Vf, 03.67.Lx 


I. INTRODUCTION 

Solitonic excitations in quantum superfluids are of sig¬ 
nificant importance. They behave like quantum mechan¬ 
ical matter waves and maintain their shape during prop¬ 
agation. As highly nonlinear and localized topological 
excitations, solitons provide a very sensitive probe of the 
fundamental coherence properties of the underlying su- 
perfluid state in which they propagate. In the ultracold 
matter of weakly interacting Bose-Einstein condensates 
(BECs), where the phase coherence is associated with 
long-range order in the one-body density matrix, solitons 
have been investigated extensively over the past decade, 
both theoretically and experimentally [If. In a series of 
ground-breaking experiments Si, dark solitons - ap¬ 
pearing as a suppression in the density profile - have 
been created via phase-imprinting. Their theoretical de¬ 
scription is provided by the nonlinear Gross-Pitaevskii 
equation p). 

In strongly interacting fermionic superfluids, solitions 
are more interesting, as the phase coherence between 
the underlying bosonic entities (i.e., Cooper pairs), char¬ 
acterized by long-range order in the two-body density 
matrix, is more subtle to understand and describe [J- 
[sj]. In addition, fermionic bound states may be induced 
near solitonic defects i, analogous to the famous An¬ 
dreev bound states inside vortex cores [l0|. In this re¬ 
spect, the recent experimental realization of dark soli¬ 
tons in strongly-interacting atomic 6 Li Fermi gases at 
the crossover from BECs to Bardeen-Cooper-Schrieffer 
(BCS) superfluids is of great interest 0 - A fermionic 
soliton was nucleated in a controlled way by using the 
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phase imprinting method in a cigar-shaped Fermi cloud 
and was observed as a reduced density slit running 
through the middle of the cloud. The soliton exhibited 
the expected oscillation when it moved from one end of 
the trap to the other. However, the rate of oscillation was 
much slower than that predicted from time-dependent 
mean-field calculations j6H9|. This discrepancy is now 
solved by a refined measurement fl2| . Due to the intrin¬ 
sic snake instability in three dimensions, the observed 
defect is actually the decay product of soliton. In the 
constrained geometry of the cloud, it is a single straight 
vortex line and therefore is better named as a solitonic 
vortex [T3 |. 

In this work, we consider the observation of dark soli¬ 
tons in an even more intriguing situation - topological 
fermionic superfluids [lj|. Our research is motivated 
by the recent opened perspective that atomic topolog¬ 
ical superfluids might be realized very soon in cold-atom 
laboratories with spin-orbit coupled Fermi gases of 40 K 
or 6 Li atoms [14! 1 1 (l | . Topological superfluids are novel 
states of matter, which have attracted considerable at¬ 
tention because of their non-trivial topological properties 
and their ability to host exotic quasi-particles known as 
Majorana fermions - particles that are their own anti¬ 
particles 00. It is therefore of great interest to ask: 
will there be any interesting features exhibited by a su¬ 
perfluid when its topological order and a solitonic defect 
come into play? 

To address this problem, we theoretically investigate 
the existence of dark solitons in a one-dimensional (ID) 
spin-orbit coupled atomic Fermi gas in a harmonic trap 
under an external Zeeman field. A dark soliton is charac¬ 
terized by a phase jump of 180° at a point node at which 
the order parameter changes sign and crosses zero. As 
the number of dark solitons can be controlled by chang¬ 
ing the number of phase jumps, we consider both cases 
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of a single soliton and multiple solitons (i.e., a soliton 
train). In the topological superfluid phase, we find that 
each dark soliton is able to host two Majorana fermions, 
well localized around nodal point of the soliton. Poten¬ 
tially, this may provide an ideal scenario to create and 
move Majorana fermions towards realistic applications, 
via the control of phase-imprinting. We propose that ex¬ 
perimentally the existence of dark solitons may be probed 
by using radio-frequency spectroscopy for the local den¬ 
sity of state or absorption imaging for the density profile. 

It should be noted that a dark soliton in one dimension 
behaves very similarly to a vortex in two dimensions. The 
latter is also a topological defect that can host a Majo¬ 
rana fermion in the vortex core and exhibit it in the local 
density of state and density profile [fH!]. A vortex lattice 
in topological superfluids has been proposed to be an 
appealing platform to perform topological quantum in¬ 
formation processing and quantum computation [ 2 p[ f 2 lj . 
In principle, a soliton train would achieve a similar goal. 
A detailed discussion of this possibility will be addressed 
in a future publication. 

We also note that recently Xu et al. investigated the 
properties of a single soliton in a ID spin-orbit coupled 
Fermi gas [22|. These authors considered a different set of 
parameters, with which the Fermi cloud enters a partial 
topological superfluid phase by increasing the external 
Zeeman field. Our results are in qualitative agreement 
with theirs when there are overlaps. 

Our paper is arranged as follows. In the next section 
(Sec. II), we briefly introduce the model Hamiltonian 
and explain how to solve it in the mean-field picture of 
the Bogoliubov-de Gennes (BdG) equations, and then 
specify the parameter space (i.e., phase diagram) for our 
numerical results. In Sec. Ill, we study the properties of 
a single soliton and show the emergence of zero-energy 
Majorana fermions when the topological regime is ap¬ 
proached. The wavefunctions of the Majorana fermions 
and their manifestations in the local density of state and 
density profile are discussed in detail. In Sec. IV, the 
cases of two and more solitons are considered. Finally, 
Sec. V is devoted to summaries and conclusions. 


II. MODEL HAMILTONIAN AND MEAN-FIELD 
THEORY 


We consider a spin-1/2 40 K Fermi gas of N atoms with 
spin-orbit coupling confined in a ID harmonic trap 
126( 1 ■ This system could be realized at Shanxi University, 
by adding a very deep 2D optical lattice (in the transverse 
y-z plane) to a spin-orbit coupled 3D Fermi gas formed by 
two counter-propagating Raman laser beams (along the 
^-direction) (14]. It can be described by a single-channel 
model Hamiltonian H = f dx[M o + J4% n t], where 


— 


— 


■^\e l2kRX ^ i + H.c. 


(1) 


is the single-particle part and 

Mnt = 9iD^\ (x) E (x) HU (x) \F t (x) ( 2 ) 

is the part describing the contact interaction between the 
two spin states. Here, a =j\ / is the pseudo-spin denoting 
the two hyperfine states and 'Fj. (x) is the fermionic field 
operator that creates an atom with mass m in the spin 
state a. The second term in describes a synthetic 
spin-orbit coupling, where flu and kn are the Rabi fre¬ 
quency and the wavevector of the laser beams, respec¬ 
tively. Due to the counter-propagating configuration, 
the momentum transferred during the two-photon Ra¬ 
man process is 2 hkn- The Hamiltonian 


Hs = 


H 2 d 2 
2m dx 2 


H + V T (x), 


(3) 


with /i being the chemical potential and Vr (x) = 
mu> 2 x 2 /2 the ID harmonic trapping potential with an 
oscillation frequency u>. The motion of atoms in the y-z 
plane is frozen due to the tight 2D optical lattice (i.e., 
the transverse trapping frequency w_l ~ Nui to). In 
this quasi-ID geometry, it is known that the low-energy 
scattering properties of atoms can be well described us¬ 
ing a contact potential giD&{x), where the ID effective 
coupling constant giD < 0 can be expressed through the 
3D scattering length 0,30 [Ij], 


2 ^ 2 «3 D _ 1 _ 

9 1D ma 2 j_ (1 — Aci3d/cl±) ’ 

where a± = y/h/(moj±) and the constant A = 
— <((l/2)/\/2 ~ 1.0326. The interatomic interaction of 
our trapped system can be conveniently parameterized 
by a dimensionless interaction parameter |28H30| 


mg id 
7 “ H 2 n 0 ’ 


(5) 


which is basically the ratio between the mean-field in¬ 
teraction energy and the kinetic energy. Here n 0 = 
(2/ir)y/Nmuj/h is the total atomic density of a non¬ 
interacting Fermi gas at the trap center within the local 
density approximation. Experimentally, in the vicinity of 
Feshbach resonances (i.e., CI3D — >• ± 00 ), the typical value 
of the interaction parameter 7 is about 5 [311. 

To solve the single-channel model Hamiltonian H, it is 
useful to apply a local gauge transformation [23j, [2J, [2(| , 


T t (x) = e + ^-L[V> t (*)-%(*)], (6) 

'M^) = e-' kRX -j=[ifr(x)+iipi(x)\. (7) 


Using the new field operators 1/7 (x) and ipi ( x ), we can 
recast the single-particle Hamiltonian into the form, 




H 0 


( x ) > (a) Ho 

_ Va (z) J ’ 

(8) 

h 2 d 2 


(9) 

~2^d^ + VT{x) 

— IX— ha z + \k x a y , 
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where we have redefined the chemical potential p, —> /i — 
H 2 ^/(2m) to absorb a constant energy shift and have 
introduced the momentum operator k x = —id/dx, the 
spin-orbit coupling constant A = h 2 ka/m and an effective 
Zeeman field h = Qr/2. Furthermore, a y and a z are 
Pauli matrices. The form of the interaction Hamiltonian 
is invariant under the local gauge transformation, i.e., 

Mnt = giD^l (x) i’l (x) 4>i ( x ) r/’f (x) ■ (10) 

The operator for the total atomic density h(x) = 
H (x) Va (x) = J2*Tpi(x)ip'T(x) is also invariant. 
However, the form of the density operator for each spin 
component changes (24| . 


A. Bogoliubov-de Gennes equations 


We solve the single-channel model Hamiltonian within 
the standard mean-field framework. By introducing an 
order parameter A ( x ) = —g id (ip^ (x) V’t ( x )), the inter¬ 
action Hamiltonian can be approximated as, 




A (x) r/'l (x) 4>l (x) + H.c. 


|AQr)| 2 

giD 


( 11 ) 


It is then convenient to use a Nambu spinor ip(x) = 
[ijjf (x ), i/.’i (x ), (x ), (a;)] T and rewrite the model 

Hamiltonian in a compact form, 


H m f = \ [ dx-ip^'HBdG'tp + Tr'Hs- f dx ^^ , (12) 
2 J J giD 

where 


' Hs — h —Xd/dx 0 — A(x) 

_ Xd/dx Hs + h A(x) 0 

riBdG = 0 A* (a:) -Hs + h Xd/dx 

—A*(x) 0 —Xd/dx —Hs — h 

. , (13) 

and the term Tr"Hs comes from the anti-commutativity of 
Fermi field operators. The mean-field model Hamiltonian 
can then be diagonalized as 


T-iBdG^r, ( x) = E v $ v (x) , (14) 

where ® v (x) = [u^ v (x) , (x), v tv (x), v± v ( x)] T and 

E v are respectively the wave-function and the energy 
of Bogoliubov quasiparticles, indexed by an integer sub¬ 
script r] = 1,2,3, The BdG Hamiltonian Eq. (fl3l) 
includes the order parameter A (x) that should be deter¬ 
mined self-consistently: 

A (z) = \ U tn V ir,f ( E v) + Hv v tvf (- E ri)] > 

(15) 

where / (E) = l/[e E ^ kBT + 1] is the Fermi distribu¬ 
tion function at temperature T. The chemical poten¬ 
tial p can be determined using the number equation, 


N = f dxn(x), where the total atomic density is given 
by 


n ( x ) — 2 f (E^) + \v ar] \ f ( E v ) 

err) 


(16) 


It is worth noting that the use of Nambu spinors double 
the Hilbert space of the system. As a consequence, we 
always have an intrinsic particle-hole symmetry in the 
Bogoliubov solutions. That is, for any “particle” solution 
with the wave-function (x) = [u^, up,, v-\- v , and 
energy E^ > 0, there is another partner “hole” solution 

with (x) = [v^,vl v ,u^ v ,ul v J T and E^ = -E^ p) < 
0. These two solutions actually correspond to the same 
physical quantum state. To avoid double counting, an 
extra factor of 1/2 appears in the expression for the order 
parameter Eq. m and the total atomic density Eq. 

m- 

The Bogoliubov equation Eq. m can be solved itera¬ 
tively with Eqs. (1151) and (fl6l) by using a basis expansion 
method, together with a hybrid strategy that takes care 
of the high-lying energy states [H, [25|, |29|, [30]. Once 
we have the solution, we calculate straightforwardly the 
local density of states, 

p (x,u) = \ y~[ [Jm CT7) | 2 5 (u — E v ) + \v av \ 2 5 (uj + E v ) . 

arj 

( 17 ) 

We note that in low dimensions phase fluctuations 
are generally enhanced due to the reduced dimenional- 
ity. As a result, in one dimension a true long-range or¬ 
der (i.e., characterized by an order parameter) is ruled 
out by the well-known MerminWVagner-Hohenberg the¬ 
orem js[i, [33}. Nevertheless, at zero temperature the 
ID pair correlation shows a slow power-law decay, much 
slower than an exponential decay expected for a normal 
state, as predicted by using bosonization or exact Bethe 
ansatz approach |3J|. Thus, we anticipate that the as¬ 
sumption of an order parameter and the use of mean- 
field BdG equations may still provide a useful approx¬ 
imation. Indeed, without spin-orbit coupling, we have 
checked that the mean-field approach leads to reason¬ 
ably accurate equation of state for a weak-coupling uni¬ 
form Fermi gas, compared with the exact Gaudin-Yang 
solution [ 29 )]. In harmonic traps, the mean-field theory 
also predicts very similar density profiles as the Bethe 
ansatz approach (within the local density approximation) 
(29[]. Therefore, in the presence of spin-orbit coupling, we 
assume that our mean-field treatment may qualitatively 
capture the essential physics of fermionic solitons. 


B. Solitonic order parameter 

A stationary dark soliton is characterized by a 7r-phase 
jump. In this case, the order parameter may be chosen 
to be real and a stationary soliton at the point node X\ 
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(* = 1,2,3, - - ■) is given by, 


A (a;) 


| A (a:)| exp 


J7T> 0 (X — Xi) 


(19) 


In numerical calculations, we calculate the order param¬ 
eter through Eq. m and then update it by phase¬ 
imprinting the 7r-phase jumps with the use of Eq. m 
or Eq. OH). The solitonic order parameter is obtained 
self-consistently, after a number of iterations up to con¬ 
vergence. 


C. Phase diagram without solitons 



0.0 0.2 0.4 0.6 0.S i.O 1.2 1.4 

h/E 
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FIG. 1: (color online) Zero-temperature phase diagram with¬ 
out solitons at 7 = 2.2 and A = 1.5Ef/ kF, determined from 
the behavior of the two lowest-energy particle solutions (with 
E ,, > 0) in the Bogoliubov quasiparticle spectrum (shown by 
the solid line and empty squares, respectively). The topologi¬ 
cal phase transition occurs at about h ~ 0.59i?F. As a result, 
two zero-energy Majorana fermions emerge at the trap edges, 
as shown in the inset, plotting the local density of states for 
the case h = 0.8-Ef- They are highlighted by two grey circles. 


may take the following form for the order parameter: 

A (x) = |A (a;)| exp [i7T0 (x — aq)], (18) 

where 0(a;) is the Heaviside step function. Similarly, the 
order parameter of a soliton train with point nodes at Xi 


The solution of the BdG equations without solitons 
has been discussed in detail in earlier works [23|, [ 2 ^, j 2 f| . 
In Fig. m we present the typical phase diagram at 
zero temperature for a spin-orbit coupled Fermi gas with 
TV = 60 atoms at the interaction strength 7 = 2.2 
and spin-orbit coupling strength A = 1.5 Ep/kF- Here, 
we use the Thomas-Fermi energy and wavevector of a 
non-interacting trapped Fermi gas, Ef = ( N/2)hut and 
= y/2 itiEf, as the units of energy and wavevector, re¬ 
spectively. The units of length are given by the Thomas- 
Fermi radius, Xf = 7 / Nh/ (mu). Throughout this work, 
we focus on the case of zero temperature and shall use 
the same total number of atoms (N = 60), interaction 
parameter (7 = 2 . 2 ) and spin-orbit coupling strength 
(A = 1.5 Ep/kF), unless otherwise specified. For the 
basis expansion, we have considered 540 single-particle 
eigenstates of the harmonic oscillator Hs and have taken 
a cut-off energy E c = 15 Ep = 450 Hu, above which we 
use the local density approximation [2J, [25| . 

Fig. H] shows the energies of the two lowest-energy 
particle solutions with E v > 0 , plotted respectively by 
using the solid line and empty squares. A topological 
phase transition occurs at a critical effective Zeeman field 
h c ~ 0.59 E F j35l - [37| . as revealed by minjl?,,} (solid line). 
At a small Zeeman field h < h c , the Fermi gas is a 
standard BCS superfluid, with a fully gapped quasipar¬ 
ticle energy spectrum (i.e., min{£' J) } > 0). Above the 
threshold, h > h c , the quasiparticle energy spectrum is 
again gapped in the bulk , as seen from the empty squares. 
However, gapless edge excitations - Majorana fermions - 
emerge at the trap edges. This is particularly evident 
when we plot the local density of states in the inset. Two 
(nearly) zero-energy Majorana fermions - well localized 
at the trap edges as highlighted by the two grey circles - 
are clearly visible. 

It is useful to note that Majorana fermions may acquire 
an exponentially small energy due to the finite size of the 
harmonic trap [38|. The typical spatial extension of the 
localized Majorana wave-function is at the order of 
the coherence length £ c = Hvf/ A, where vf and A are 
the unperturbed local Fermi velocity and pairing gap at 
the trap edge (H). For the two Majorana fermions shown 
in the inset of Fig. [T] we estimate that ~ O.lzj? (see 
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FIG. 2: (color online) (a) The solitonic pairing gap A(x) at 
zero Zeeman field h = 0 and at different interaction strengths 
7 and spin-orbit coupling strengths A, when a 7r-phase jump 
is imprinted at xi = 0. (b) The solitonic A(x) at nonzero 
Zeeman fields [h/E F = 0.4 (thick solid line), 0.56 (crosses), 
0.58 (daggers) or 0.8 (dashed line)] and at 7 = 2.2 and A = 
1.5Ef/ k,F ■ For the case of h = QAEf, we also consider a 
soliton at *1 = — 0.5:ef and plot the result with a thin solid 
line. 

also Fig. 0 . Therefore, the exponentially small overlap 
between two Majorana fermion wave-functions and 
i.e., O = ~ exp(—L/£ m ), where L ~ 2x F 

is the distance between two Majorana fermions, leads 
to an exponentially small energy (splitting) of Majorana 
fermions, 

it (20) 

This is consistent with our numerical finding that the 
energy of Majorana fermions E ~ 10” 1 o £’f- 

III. SINGLE SOLITON 

Here we consider the behavior of a single dark soliton. 
The case of multiple dark solitons will be discussed in the 
next section. 


A. Order parameters 

In Fig. [51 we report the pairing gap profile in the pres¬ 
ence of a single dark soliton, at zero Zeeman field (a) 
or finite Zeeman fields (b). The order parameter crosses 
zero at the position of the soliton x\ and hence creates 
a point node. At small Zeeman fields, it exhibits two 


length scales around the point node [l|: a fast oscillation 
with length scale of kp , and a slower healing with length 
scale £ c = Kv F / A. Here, v F and A are the unperturbed 
local Fermi velocity and pairing gap at the point node X \, 
respectively. The former length scale is essentially inde¬ 
pendent of the interaction parameter and the spin-orbit 
coupling strength. Thus, as in the case of a vortex in 2D 
Fermi gases [39|, we may safely identify the oscillation as 
the Friedel oscillation. For the coherence length, we find 
that £ c ~ 3kp 1 at 7 = 2.2. It increases with decreasing 
interaction parameter, as expected. 

For a spin-orbit coupled Fermi gas, it is interesting 
to see (Fig. [2jb)) that the Friedel oscillation suddenly 
ceases to exist when the Zeeman field is above a thresh¬ 
old, h c — 0.57 E f . Actually, this point corresponds to 
the topological phase transition in the presence of a sin¬ 
gle dark soliton, which we shall now discuss in greater 
detail. 


B. Andreev bound states and Majorana fermions 

In the presence of a single dark soliton, the topolog¬ 
ical phase transition point can also be determined from 
the Bogoliubov quasiparticle spectrum. In Fig. 0a), we 
plot the energies of the three lowest-energy particle so¬ 
lutions as a function of the Zeeman field, while in Figs. 
0 b) and 0 c), we show the characteristic spectrum be¬ 
fore and after the topological transition. The transition 
is associated with the closing and reopening of the energy 
gap in the bulk [H, [H, IH, H 3 , HE [36|, see, for example, 
the empty squares in Fig. 0a) for the lowest-energy bulk 
state in the particle branch. Thus, we determine that the 
transition occurs at the critical field h c — 0.57E F , which 
is slightly smaller than the threshold h c ~ 0.59E F in the 
absence of the dark soliton (see Fig. 0). 

In the topologically-trivial BCS superfluid phase, we 
find two Andreev-like bound states that reside near the 
point node of the soliton. Their existence is easy to un¬ 
derstand from the density dip of the solitonic order pa¬ 
rameter, which basically creates an effective confinement 
potential of length scale £ c for quasiparticles. As a result, 
localized states develop, with a characteristic energy of 
the order of K 2 /(mt; 2 ) = A 2 /(2e^), where e F is the local 
Fermi energy. These are reminiscent of the well-known 
Caroli-de Gennes-Matricon states in a vortex core in 2D 
Fermi gases 0 - 

At zero Zeeman field, the two Andreev bound states 
are degenerate in energy, which is nonzero and is called 
the minigap in the context of superconductors in solid 
state physics. In the weakly-interacting limit, as shown 
in the inset of Fig. 0 a), we have checked that, to a very 
good approximation, 

A 2 

^minigap ^ 0.26 — . ( 21 ) 

£f 

With increasing Zeeman field, the energies of the two An¬ 
dreev bound states initially split: one increases and the 
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FIG. 3: (color online) (a) The energies of the three lowest-energy particle solutions with Erj > 0, as a function of the Zeeman 
field in the case of a single dark soliton at xi = 0. The solid and empty circles show Andreev bound states localized near the 
soliton, while the empty squares correspond to the lowest-energy particle state in the bulk. The vertical dashed line separates 
the BCS and topological superfluid phases. The inset shows the minigap at h = 0 and at different interaction parameters 
7 and spin-orbit coupling strengths [XIcf/Ef = 1.5 (circles), 1.0 (squares) and 0.0 (daggers)], as a function of the square of 
the local pairing gap, A 2 /e%. The dashed line is the linear fit -Emmigap = 0.26A 2 /tF- (b) and (c) The spatial distribution of 
Bogoliubov quasiparticle energy spectrum at h = OAEf and h = 0.8Ef, respectively. The Andreev bound states near the 
soliton are highlighted by empty and solid circles. Here, we approximately characterize the location of a quasiparticle by using 
its wave-function: yj( x 2 ) = {f dxx 2 X] ct [ w ct i x ) + (a;)]} 1 / 2 , if the wave-function is well-localized at a certain point. Otherwise, 
a/( a: 2 ) characterizes the width of the wave-function. 


other decreases. However, when h > 0.2 Ep, both ener¬ 
gies gradually decrease to zero as h nears the topological 
transition. A typical energy spectrum before the tran¬ 
sition at h = OAEf is shown in Fig. [31b), plotted as a 
function of the approximate location of each quasiparticle 
state, 


VW) = 


1 1/2 


dx x 2 V (u 2 (x) 




0 )) 


( 22 ) 


The two Andreev bound states, represented by the empty 
and solid circles [40], are clearly seen near the point node 
of the soliton at X\ = 0. The highly-localized wavefunc- 
tions of the lowest-energy bound state, corresponding to 
the solid circle, are plotted in Fig. [4|a). 

Approaching the critical Zeeman field h c ~ 0.57 Ep, 
the energies of the two Andreev bound states tend to 
zero. The energies of some bulk states also vanish. The 
interference of these nearly zero-energy modes leads to a 
huge reconstruction of the quasiparticle spectrum across 
the topological transition. Immediately after the transi¬ 
tion, with the reopening of the energy gap in the bulk, 
we observe that one of the Andreev bound states merges 
with bulk states and therefore can no longer be identified 
as a localized state. At the same time, two zero-energy 


edge states appear at the trap edges. As a result, we 
find four Majorana fermions when the system is in the 
topological phase: two at the edges and the other two 
near the soliton [4p| |. Physically, the number of Majo¬ 
rana fermions may be understood from the fact that a 
single dark soliton effectively splits the Fermi gas into 
two, each of which could host two Majorana fermions. 
In the case of multiple dark solitons, we therefore antic¬ 
ipate that the total number of Majorana fermions would 
be 2n + 2, where n is the number of solitons. This ex¬ 
pectation is confirmed by our numerical calculations. 

In Figs. Htb) and HJc), we examine the wave func¬ 
tions of the four Majorana fermions. Since all four states 
are degenerate with zero energy, they may be mutually- 
superposed. The interference leads to very similar wave- 
functions, so in the figure, only one of the four is plot¬ 
ted. The bond and anti-bond superpositions of the well- 
localized Majorana wavefunctions, presumably one at 
the point node of the soliton and the other two at the 
edges, are fairly clear (see the next paragraph for more 
discussions) [22|. Each of these localized wavefunctions 
satisfies the symmetry of u a ( x ) = ±i/* ( x ), which is pre¬ 
cisely the required symmetry for Majorana fermions. As 
a result, the approximate location yj(x 2 ) of the four Ma¬ 
jorana fermions is exactly the same (i.e., ill-defined), as 
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FIG. 4: (color online) The wavefunctions of the lowest-energy 
Andreev bound states in the BCS superfluid phase at h = 
QAEf (a) or in the topologically-nontrivial superfluid phase 
at h = 0.8 Ef (b and c). In the latter case, the wavefunctions 
strongly interfere with that of Majorana fermions at the trap 
edges. The superposition persists when we place the dark 
soliton in a less symmetric position x\ = — O.Ixf- 

shown in Fig[3](c). This superposition is very robust with 
respect to the position of the dark soliton, as can be seen 
in Fig. He). When we displace the dark soliton away 
from the origin to the left ( x\ = — O.la;^), the superpo¬ 
sition remains, but the relative weights of the Majorana 
wavefunctions localized at the two edges are no longer 
equal. 

It is important to note that, despite the superpo¬ 
sition in wavefunctions, the energies of all the Majo¬ 
rana fermions are nearly zero (i.e., about 10~ 1O -Ef in 
our numerical calculations) [22j]. This fact is consistent 
with the earlier observation that the only one Majo¬ 
rana fermion wave-function at the point node of the soli- 
ton is superposed with the two edge Majorana fermions, 
so that the energy splitting is exponentially small as 
in Eq. ®. Otherwise, if the two solitonic Majo¬ 
rana fermions near the point node do interfere with each 
other, we would rather anticipate a sizable energy split¬ 
ting E/Ep ~ e~ L ^ M ~ 1, since now the distance L ~ 0. 


C. Experimental detection of Majorana fermions 

We now turn to consider the experimental observa¬ 
tion of the two additional Majorana fermions localized 
at the point node of the dark soliton. A direct and con¬ 
venient way is to use spatially-resolved radio-frequency 
(rf) spectroscopy, which acts as a cold-atom analog of 
scanning tunneling microscopy (STM) and measures the 
local density of states j4j, |4j|. In Fig. H we show the 
local density of states p(x, w) before and after the topo- 
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FIG. 5: (color online) Local density of states of a BCS su¬ 
perfluid at h = QAEf (a) and of a topologically-nontrivial 
superfluid at h = 0.8 Ef (b), when a single dark soliton is cre¬ 
ated at x\ = 0. In the topological phase, the two Majorana 
fermions near the point node of the dark soliton (at xi = 0) 
are not distinguishable because of the superposition of the 
wavefunctions. The color map indicates the magnitude of the 
local density of states in units of tif/Ef- 


logical transition. In the BCS superfluid phase (a), the 
solitonic Andreev bound states can be easily identified. 
In the topologically-nontrivial phase (b), we can clearly 
see the two zero-energy Majorana fermions residing at 
the two trap edges. In addition, there is a zero-energy 
response near the origin, arising from the two Majorana 
fermions localized near the dark soliton. However, due 
to their superposition (i.e., overlapping wave functions), 
they are not individually resolvable. 

Alternatively, the existence of the soliton-induced Ma¬ 
jorana fermions may be indirectly deduced from the to¬ 
tal density profile, which could be measured via in-situ 
or time-of-flight absorption imaging. In Fig. [Gl we 
present the density profile at two Zeeman field-strengths. 
Before the topological transition at h = OAEp (solid 
line), there is an apparent oscillation with spatial period 
~ kp , in accord with the Friedel oscillation in the soli¬ 
tonic order parameter (see Fig. H. As we approach the 
topological transition point at the critical field strength 
h c ~ 0.57Ep, the amplitude of the density oscillation 
An = n max — n m i n , where n max and n m i n are respectively 
the maximum and minimum densities near the soliton, 
rapidly decreases and vanishes precisely at the transition 
(see the top right inset of Fig. H. After the topological 
transition (see, for example, the dashed line at h = 0.8 Ep 
in the main figure), the density profile becomes flat and 
the peak density at the trap center is essentially inde¬ 
pendent of the Zeeman field. The disappearance of the 
density oscillation is associated with the formation of the 
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FIG. 6: (color online) The total density profile in the presence 
of a single dark soliton at x\ = 0. The solid and dashed lines 
show the results at h = OAEf and h = 0.8 Ef, respectively. 
The inset at the bottom is an enlarged view around the origin 
x = 0, with additional results near the topological transition 
at h = 0.56-Ef (crosses) and at h = 0.58 Ef (daggers). The 
inset at the top right top shows the oscillation amplitude An 
of the density profile near the origin, An = Umax — n m i n , where 
Umax and n m i n are the maximum and minimum densities, re¬ 
spectively. 


solitonic Majorana fermion modes, whose occupation sig¬ 
nificantly contributes to the total density because of the 
large amplitude of their localized wavefunctions (see Fig. 
[3Jb)). This is very similar to what happens in the vor¬ 
tex core of a topological superfluid |19j. There the core 
density is also greatly affected by the formation and oc¬ 
cupation of the vortex-core Majorana fermion modes. 

In connection with current experiments, we may con¬ 
sider a spin-orbit coupled 40 K Fermi gas in the presence 
of a tight 2D optical lattice, with an axial trapping fre¬ 
quency u> = 2n x 116 Hz 0- For the typical num¬ 
ber of total atoms N = 60 in each tube of a ID Fermi 
cloud, the Fermi temperature Tp is about 200 nK. We 
may take kp — 0.75 kp and then the recoil energy is 
Ek ~ 0.56Ep. The topological phase transition at low 
temperatures (i.e., T < 0.1 Tp ~ 20 nK) takes place at 
the critical Zeeman field h c ~ 0.6 Ep ~ Ep, which cor¬ 
responds to a critical Rabi frequency flp ~ 2 Ep. The 
size of the Fermi cloud is about xf ~ 12 fim. To resolve 
the zero-energy Majorana modes by rf spectroscopy, we 
may require the frequency/energy resolution to be better 
than 27t x 100 Hz. On the other hand, for the in-situ 
density profile, the spatial resolution needed to measure 
the amplitude of the density oscillations is about 0.5/nm. 
Thus, it seems more practical to use absorption imaging 
after some time-of-flight, if we assume that the structure 
of density oscillations survives for a short expansion time. 



x/x ¥ 



FIG. 7: (color online) The solitonic order parameter (a) and 
the total density distribution (b) in the presence of two dark 
solitons (placed at xi = — 0.5:cf and *2 = +0.5:ef), at two 
Zeeman fields: h = OAEf (solid line) and h = 0.8 Ef (dashed 
line). 


IV. MULTIPLE SOLITONS 

We now consider a soliton train. Without loss of gener¬ 
ality, let us focus on the case of two dark solitons. Other 
cases with three or more dark solitons were also examined 
and so far we have not found additional, qualitatively dif¬ 
ferent results. 

In Fig. [71 we plot the solitonic order parameter and 
the total density distribution, when the two solitons are 
placed at X\ = —0.5a;F and X 2 = +0.5:z;,f- Qualitative 
features such as the Friedel oscillations in the order pa¬ 
rameter and density profile near the point nodes of the 
solitons in the BCS superfluid (h = 0 AEp, solid line), 
can be easily understood analogously to the case of a 
single dark soliton, as discussed earlier (cf. Figs. [2] and 
©. 

In Fig. [51 we present the wavefunctions of three Majo¬ 
rana fermions in the topological phase (a, b and c) and 
their manifestation in the spatially-resolved rf spectrum 
(d). In total, there are six Majorana fermions, localized 
pair-wise at the trap edges and at the point nodes of the 
solitons. Only three of them are shown, one out of each 
pair, owing to the particle-hole symmetry. Unlike in the 
case of a single soliton, the two Majorana edge inodes do 
not interfere with the solitonic Majorana fermions and 
hence have essentially zero energy (i.e., E ~ 10 -1 o .Ef) 
due to the exponentially small overlap in their wavefunc- 
tions. In contrast, the overlap of the two solitonic Ma¬ 
jorana fermion wavefunctions - originating from the two 
solitons - is significant. This leads to a nonzero energy 
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FIG. 8: (color online) (a, b and c) The wave-functions of three Majorana fermions of a topological superfluid at h — 0.8 Ef, in 
the presence of two dark solitons at xi = —0.5 xf and X2 = +0.5:tf. The energy of each Majorana fermion is indicated, (d) 
The corresponding local density of state within the bulk energy gap. The color indicates the magnitude of the local density of 
states, see, for example, Fig. 0 


for each solitonic Majorana fermion, which, following Eq. 
(1201), is of the order of 

<23) 

Here, in the estimation, we have used the distance L ~ 
xf and a bit larger localization length scale ~ 0.2 xf 
for the two solitonic Majorana wave-functions shown in 

Figs. [Sib) and [5]( c )- 

V. CONCLUSIONS 

To summarize, we have investigated the behavior of 
dark solitons in a one-dimensional topological super¬ 
fluid, in the context of ultracold atomic Fermi gases 
with _spin-orbit coupling and an external Zeeman field 
[3 HI. We have predicted that each dark soliton 
can host two Majorana fermions localized at its point 
node, which are detectable through the techniques of 
spatially-resolved radio-frequency spectroscopy or ab¬ 


sorption imaging. Therefore, the well-known technique 
of creating dark solitons via phase imprinting also allows 
one to create solitonic Majorana fermions. These Majo¬ 
rana fermions can then find realistic applications in, e.g., 
topological quantum information processing and quan¬ 
tum computation. This scheme is very similar to the idea 
of using a vortex lattice in a two-dimensional topological 
superfluid, where the Majorana fermions at the vortex 
cores are used as qubits [2l|. In current cold-atom ex¬ 
periments, it would be much easier to engineer a soliton 
train than to create a vortex lattice. For the purpose 
of exchanging solitons at different positions to demon¬ 
strate the non-Abelian statistics of Majorana fermions, 
in future studies it would be interesting to understand 
traveling grey solitons characterized by a complex order 
parameter and nonzero velocity [Q-Q- 
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